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Gene regulatory circuits must contend with intrinsic noise that arises due to finite numbers of 
proteins. While some circuits act to reduce this noise, others appear to exploit it. A striking example 
is the competence circuit in Bacillus subtilis, which exhibits much larger noise in the duration of its 
competence events than a synthetically constructed analog that performs the same function. Here, 
using stochastic modeling and fluorescence microscopy, we show that this larger noise allows cells to 
exit terminal phenotypic states, which expands the range of stress levels to which cells are responsive 
and leads to phenotypic heterogeneity at the population level. This is an important example of how 
noise confers a functional benefit in a genetic decision-making circuit. 


I. AUTHOR SUMMARY 

Fluctuations, or “noise”, in the response of a system is 
usually thought to be harmful. However, it is becoming 
increasingly clear that in single-celled organisms, noise 
can sometimes help cells survive. This is because noise 
can enhance the diversity of responses within a cell pop¬ 
ulation. In this study, we identify a novel benefit of noise 
in the competence response of a population of Bacillus 
subtilis bacteria, where competence is the ability of bac¬ 
teria to take in DNA from their environment when under 
stress. We use computational modeling and experiments 
to show that noise increases the range of stress levels for 
which these bacteria exhibit a highly dynamic response, 
meaning that they are neither unresponsive, nor perma¬ 
nently in the competent state. Since a dynamic response 
is thought to be optimal for survival, this study suggests 
that noise is exploited to increase the fitness of the bac¬ 
terial population. 


II. INTRODUCTION 

Snapshots of bacterial populations often reveal large 
phenotypic heterogeneity in the gene expression states of 
its composite individuals. Such phenotypic heterogene¬ 
ity in a clonal population of bacterial cells in a single 
environment has significant consequences for how well 
the organisms can adapt and survive. On the one hand, 
a population with little or no heterogeneity may allow 
for all cells to take advantage of certain optimal condi- 
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tions to which the population is exposed. In this case, 
heterogeneity is suboptimal and therefore detrimental to 
fitness. On the other hand, numerous recent studies have 
shown that heterogeneous populations allow for cells to 
account for uncertainty in future environmental condi¬ 
tions [1-6]. In this case, heterogeneity is benehcial to 
fitness. A straightforward way to maintain high pheno¬ 
typic heterogeneity is for each cell to exhibit a dynamic 
response. This allows each cell in its turn to transi¬ 
tion among the various states of the population, e.g. via 
switching, pulsing, or oscillatory dynamics. The hetero¬ 
geneity is intrinsically encoded in each cell, and is often 
enhanced by, or even entirely due to, stochasticity, or 
“noise”, at the molecular level [2-4, 6]. 

The ability of molecular noise to cause stochastic phe¬ 
notype changes has been demonstrated in a number of bi¬ 
ological systems. In the context of enzymes, several stud¬ 
ies have explored how intrinsic noise due to low numbers 
of molecules, or even a single molecule, can have dramatic 
effects through the amplified actions of a few enzymes 
[7, 8]. Moreover, studies of bacterial operons, including 
in the context of bacterial persistence, have suggested 
that stochasticity could be encoded in the interactions 
between genes in a genetic regulatory network by ensur¬ 
ing that certain operon states are exposed to low num¬ 
bers of molecules [9-11]. Recently, a theoretical study 
has demonstrated the conditions for when deterministic 
approaches to modeling genetic circuit dynamics break 
down, due to amplified effects of rare events caused by a 
small number of regulators [12]. Together, these works 
suggest that phenotypic heterogeneity could be rooted in 
low-molecule-number noise, and that this noise could in 
turn be encoded in the architecture of genetic regulatory 
networks. 

The competence response of the gram-positive bac¬ 
terium Bacillus subtilis provides a striking example of 
dynamically maintained phenotypic heterogeneity. Un- 
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FIG. 1: Schematic illustrating phenotypic heterogeneity and the effects of noise. (A) When all cells in a popnlation 
exhibit either no response (left) or a high response (right), then the population is homogenous. In contrast, if individual cells 
exhibit a dynamic response (middle), this leads to a heterogenous population, with a fraction / of cells in the responsive state 
at any given time. (B) Intrinsic noise affects the dynamics of the response. For the B. subtilis competence response, we find in 
this study that noise expands the viable response range: the range of stress levels over which / remains neither 0 nor 1. 


der stress, B. subtilis undergoes a natural and transient 
differentiation event, termed competence, that allows the 
organism to incorporate exogenous genes into its genome. 
Previous studies have shown that entry into the com¬ 
petent state is controlled by a genetic circuit that that 
can be tuned to one of three dynamical regimes [13]: an 
excitable regime at low stress levels, where cells rarely 
and transiently enter the competent state; an oscillatory 
regime at intermediate stress, where cells oscillate in and 
out of the competent state; and a mono-stable regime at 
high stress, where cells remain in the competent state. 
Importantly, oscillatory (and repeatably excitable) dy¬ 
namics lead to phenotypic heterogeneity, since cells are 
dynamically transitioning in and out of the competent 
state (see Fig. lA). This heterogeneity is especially im¬ 
portant to the survival of B. subtilis: if no cells respond, 
competence is not exploited, and the population may suc¬ 
cumb to the stress. On the other hand, if all cells are 
permanently in the competent state, this can also be fa¬ 
tal to the population, since competence has been shown 
to reduce the cell growth rate and prevent cell division 
due to the inhibition of FtsZ [14, 15]. Therefore, main¬ 
taining a dynamic competence response, and therefore 
a heterogenous population, is thought to be crucial to 
survival under stress. 

The effects of noise on the dynamics of the compe¬ 
tence response are only partially understood. Previous 
work has shown that noise can trigger excitations into the 
competent state when the circuit is tuned to the excitable 
regime [15]. Later work showed further that these excita¬ 
tions have a high variability in their duration, and that 
this variability is directly linked to the architecture of 


the competence circuit [16]. In particular, this work em¬ 
ployed an analogous synthetic excitable circuit, termed 
SynEx, to provide evidence that the duration variability 
is due to intrinsic noise from low molecule numbers in 
the native circuit. However, the ability of this intrinsic 
noise to trigger sustained or repeatable excitations has 
not yet been quantified. Moreover, the generic effects of 
intrinsic noise on the three dynamic regimes, and how 
these effects translate to the physiological function of B. 
subtilis at the population level, are unknown. 

Here, using stochastic modeling and quantitative flu¬ 
orescence microscopy, we study the effects of intrinsic 
noise on the competence dynamics and the ensuing pop¬ 
ulation heterogeneity of B. subtilis. We uncover a novel 
effect of noise that goes beyond architecture-dependent 
stochastic effects in a single cell. Specifically, we find 
that at both low and high stress levels, noise prevents 
cells from becoming unresponsive or indefinitely respon¬ 
sive to the stress, and instead allows cells to respond dy¬ 
namically. These effects expand the range of stress levels 
over which the population of cells maintains a heteroge¬ 
neous response distribution, which is critical to the pop¬ 
ulation viability (see Fig. IB). The use of efficient numer¬ 
ical methods and stochastic simulation at several levels 
of model complexity allows us to elucidate the mecha¬ 
nisms behind these effects. A central prediction from our 
modeling is that these effects are rooted in noise arising 
from low numbers of molecules. We verify this prediction 
using quantitative fluorescence microscopy by comparing 
the population response of native B. subtilis with that of 
synthetic mutants harboring the less-noisy SynEx circuit. 
Taken together, these results constitute a fundamental 
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FIG. 2: Architectures and model parameters of the native and SynEx circuits. The top row summarizes the 
regulatory interactions, while the bottom row depicts the model details. (A) In the native circuit, ComK is produced with the 
induction rate and activates its own expression with Hill function parameters /3fc, kk, and h. ComS is expressed at the basal 
rate Os and is repressed by ComK with Hill function parameters /3s, ks, and p. ComK and ComS are degraded at rates A*, 
and As, respectively, and, additionally, both compete for binding to the degradation enzyme MecA. MecA degrades ComK and 
ComS with maximal rates and ^s, respectively, and with Michaelis-Menten constants Tk and Fs, respectively. (B) In the 
SynEx circuit, ComK is produced with the induction rate ak and activates its own expression with Hill function parameters 
/3fe, kk, and h. MecA is expressed at the basal rate am and is activated by ComK with Hill function parameters /3m, km, and 
p. ComK and MecA are degraded at rates A*, and Am, respectively, and MecA enzymatically degrades ComK with rate S. 


example of how noise can increase the functionality of a 
phenotypic response. 


III. RESULTS 

Entry of B. subtilis cells into the competent state oc¬ 
curs at high expression levels of the ComK protein. This 
protein activates a set of downstream genes allowing for 
the uptake of DNA [15]. ComK is typically expressed at a 
basal level, and stress in the environment alters the level 
of expression. In our genetic circuit design, as described 
below, increasing the stress level is mimicked by induc¬ 
ing comK expression using an increasing amount of a 
lactose analogue. Isopropyl /3-D-l-thiogalactopyranoside 
(IPTG), in the environment. 

In the native competence circuit, ComK activates its 
own expression, and represses the expression of another 
protein, ComS. ComS and ComK compete to be rapidly 
degraded by the MecA protein complex [15] (see Fig. 2A, 
bottom). Therefore, high concentrations of ComS hinder 
the degradation of ComK, effectively providing positive 
feedback to ComK by allowing ComK levels to build up. 
These interactions are summarized in Fig. 2A (top). 

In the SynEx circuit, as described in [16], the repres¬ 
sion of ComS by ComK is removed by gene knockout. 


Then, the expression of MecA is placed under the control 
of ComK. This causes ComK to activate MecA, which in 
turn represses ComK via active protein degradation (see 
Fig. 2B, bottom). These interactions are summarized in 
Fig. 2B (top). Note that in the native circuit, ComK 
represses its own activator (ComS), while in the SynEx 
circuit, ComK activates its own repressor (MecA). 


Both the native and SynEx circuits have architectures 
characteristic of molecular oscillators. Therefore we ex¬ 
pect both circuits to allow for a dynamic response of each 
individual in a population. However, the main difference 
is that in the native circuit, when ComK levels are high, 
ComS levels are low, which leads to large amounts of 
intrinsic noise. In contrast, in the SynEx circuit, when 
ComK levels are high, MecA levels are also high, corre¬ 
sponding to less intrinsic noise. Previous work showed 
that this difference in architecture causes the native cir¬ 
cuit to display a broad range of competence durations, 
whereas the SynEx circuit displays a relatively narrow 
range of competence durations [16]. However, the effects 
of noise and architecture on the ranges of dynamic re¬ 
sponse and the ensuing population heterogeneity in these 
systems remained unknown. 
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FIG. 3: Stochastic modeling of B. subtilis competence. (A) The deterministic model of each circuit exhibits three dynamic 
regimes (excitable, oscillatory, and mono-stable), depending on the ComK induction rate Ok, which models stress level. (B) 
The stochastic model reveals the ensuing distribution of response levels in each of the dynamic regimes. The fraction of the 
distribution in the responsive state / (determined by the inflection points, see Sec. V) is shaded. (C) Whereas the deterministic 
model exhibits sharp transitions between the dynamic regimes (dashed lines), the stochastic model exhibits a continuous 
dependence of / on induction rate. We see that for both circuits, stochasticity extends the viable response range (0 < / < 1) 
beyond the transitions predicted by the deterministic model, in both directions. Parameters are as in [16] and are given in 
Appendix A. In A and B, from left to right, the values of the control parameter are a*, = {0.072,1.15, 36}/hour (native) and 
ak = {0.036,1.8,36}/hour (SynEx). 


A. Noise expands the response range 

To elucidate the effects of noise in each of the na¬ 
tive and SynEx circuits (Fig. 2), we develop a stochas¬ 
tic model of each circuit, which includes noise, and then 
compare each to its deterministic analog, which does not 
include noise. As described in Sec. V, Materials and 
Methods, we develop the stochastic models at several 
levels of complexity to investigate the robustness of our 
findings to our modeling assumptions, and we solve each 
model using a combination of efficient numerical solution 
and stochastic simulation. We first describe the behavior 
of the deterministic models. As shown in Fig. 3A, a stan¬ 
dard linear stability analysis of the deterministic model 
for each circuit reveals three dynamical regimes, depend¬ 
ing on the value of the control parameter, the ComK 
induction rate ak- At low induction, each circuit is ex¬ 
citable, resulting in a transient differentiation event into 
and out of the competent (high-ComK) state. At inter¬ 
mediate induction, each circuit is oscillatory, periodically 


entering and exiting the competent state. At high induc¬ 
tion, each circuit is mono-stable, staying in the compe¬ 
tent state indefinitely. These three dynamical regimes 
have been confirmed in experimental studies of the na¬ 
tive competence circuit [13]. 

We find that these deterministic dynamics are reflected 
in the stationary solutions to the minimal stochastic 
models. As shown in Fig. 3B, the three types of dynam¬ 
ics correspond to three shapes of stationary probability 
distributions of ComK levels. Excitable dynamics corre¬ 
spond to a distribution confined to low ComK molecule 
numbers, oscillatory dynamics correspond to a distribu¬ 
tion mixed between low and high molecule numbers, and 
mono-stable dynamics correspond to a distribution cen¬ 
tered at high molecule numbers. As described in Sec. V, 
we calculate the fraction / of the distribution in the high- 
molecule-number state (see the shaded regions in Fig. 
3B). Within our model, / represents the fraction of time 
a single cell spends in the competent state, or equiva¬ 
lently, the fraction of an isogenic population of cells found 
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FIG. 4: Stochastic oscillations persist outside the deterministic oscillatory regime. The deterministic oscillatory 
regime is defined by for the induction rate a^. (A) At low induction rate Ofe < ce\^\ where the deterministic 

model predicts excitable dynamics, the stochastic dynamics are oscillatory. The oscillations arise from repeated noise-induced 
excitations. (B) At high induction rate au > where the deterministic model predicts mono-stable dynamics, the stochastic 
dynamics are also oscillatory. The oscillations here arise because noise prevents damping to the mono-stable state (see the 
deterministic curves in the right panels). The effect is much stronger for the native circuit (notice that the left panel is 15 times 
outside the deterministically oscillatory regime) because, unlike in the SynEx circuit, one of the species, ComS, is at low copy 
number and therefore subject to significant intrinsic noise. 


in the competent state at a given time. Importantly, / 
is the indicator of population heterogeneity, since unre¬ 
sponsive (/ = 0) or fully competent (/ = 1) populations 
are homogeneous, while mixed populations (0 < / < 1) 
are heterogeneous. We define the range of induction rate 
Ofc for which 0 < / < 1 as the viable response range, 
since unresponsive cells (/ = 0) do not benefit from com¬ 
petence, while long-term competence (/ = 1) is known to 
have a detrimental effect on growth rate and cell division 
[14, 15]. 

In Fig. 3C, we compare the viable response range of the 
stochastic model with the boundaries between dynami¬ 
cal regimes predicted by the deterministic model. We 
see that for both the native and the SynEx circuit, the 
stochastic range extends beyond the deterministic range 
for both low and high induction rate a^. Furthermore, 
the extension at high induction rate is significantly more 
pronounced for the native circuit (roughly 20 times the 
deterministic value) than for the SynEx circuit (roughly 
3 times the deterministic value). These observations im¬ 
ply that noise expands the range of stress levels to which 
cells can respond in a dynamic way. In the next section, 
we elucidate the mechanisms behind this expansion. 


B. Noise-induced oscillations underlie the 
expansion of the response range 

Why does noise expand the viable response range at 
low induction levels? As shown in Fig. 4A, the reason 


is that noise leads to repeated excitations into the com¬ 
petent state, which prevents the system from remaining 
completely unresponsive. In a completely deterministic 
excitable system, an excitation is caused by initializing 
the system away from its stable fixed point, and it occurs 
only once. However, in a stochastic system, noise can 
cause repeated perturbations away from the stable state, 
leading to persistent additional excitations. Indeed, in 
both circuits, noise at the stable state is high, because 
the stable state corresponds to one or more species being 
expressed at very low molecule number (ComK for the 
native circuit, ComK and MecA for the SynEx circuit; see 
Fig. 4A). Since the dynamics are governed by Poissonian 
birth-death reactions, low molecule numbers correspond 
to high intrinsic noise (variance over the squared mean), 
leading to frequent and persistent excitations. This ef¬ 
fect is consistent with the noise-induced excitations seen 
for these circuits in previous work [15, 16]. Here, how¬ 
ever, we have quantified the effect of these excitations on 
the stochastic distribution, which describes the heteroge¬ 
neous population response. 

Why does noise expand the viable response range at 
high induction levels? Here the mechanism is different 
from at low induction levels. As shown in Fig. 4B, the 
reason is that noise prevents the damping of oscillations, 
which keeps the system from relaxing to the competent 
state. In the deterministic system, the mono-stable state 
is defined by a stability matrix whose eigenvalues are 
complex with negative real parts (Fig. Al). This means 
that the solution relaxes to the mono-stable state in an 
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oscillatory way, i.e. the oscillations are damped (see the 
black lines in the right panel of Fig. 4B, for example). In¬ 
trinsic noise thwarts this relaxation, continually perturb¬ 
ing the system away from the stable point, and preserv¬ 
ing a finite oscillation amplitude (see the colored lines). 
Similar effects have been observed in ecological and epi¬ 
demic models, where they are attributed to the ability of 
white noise to repeatedly excite a system at its resonant 
frequency [17]. Here we see the effect at the molecular 
level in bacteria, and we find that it occurs sufhciently 
strongly that it supports and significantly extends a het¬ 
erogeneous population response. 

At high induction levels, the expansion of the viable 
response range is more pronounced in the native circuit 
than in the SynEx circuit. This effect was demonstrated 
at the population level in Fig. 3C. It is also demonstrated 
by the dynamics in Fig. 4B: the noise-induced preven¬ 
tion of damping is clearly evident for the native circuit, 
even at the Uk value shown, which is 15 times value pre¬ 
dicted deterministically. The reason that the effect is so 
pronounced in the native circuit is that the mono-stable 
fixed point corresponds to ComS being expressed at very 
low molecule numbers, where the intrinsic noise is high 
(lower left panel). In contrast, in the SynEx circuit, the 
mono-stable state corresponds to both species begin ex¬ 
pressed at higher molecule numbers, so the intrinsic noise 
is lower. This difference, which stems ultimately from the 
difference in the architecture of the two circuits (Fig. 2), 
was found in previous work [16] to be responsible for the 
increased variability in the competence durations of the 
native circuit compared to the SynEx circuit. Here we 
demonstrate that the architecture of the native circuit 
additionally leads to an increase in the expansion of its 
viable response range, which has a clear benefit for fit¬ 
ness. 

We have tested that the effects discussed above are ro¬ 
bust, in that they persist when we relax the three simpli¬ 
fying assumptions of our minimal stochastic model (see 
Sec. V). We relax two of the assumptions by considering 
a non-adiabatic stochastic model in which the fast dy¬ 
namics of mRNA production and enzymatic degradation 
are included explicitly, and by setting the mean molecule 
numbers in the tens of thousands as opposed to tens (see 
Appendix A). We find that all noise-induced effects per¬ 
sist, namely (i) repeated excitations, (ii) the prevention 
of damping, and (iii) the enhancement of effect ii in the 
native circuit over the SynEx circuit (see Fig. A5 and Fig. 
A6). As shown in Fig. A5 and Fig. A6, we also verified 
quantitatively that effects i and ii produce sufficiently 
oscillatory dynamics that the power spectrum is peaked, 
as opposed to the non-peaked power spectrum observed 
for purely excitable or mono-stable dynamics. Interest¬ 
ingly, when we raise the molecule number, but retain the 
adiabatic assumption, we find that the effects of noise di¬ 
minish, and the stochastic model behaves like the deter¬ 
ministic model (see Fig. A7 and Fig. A8). This confirms 
that the effects we observe are rooted in the intrinsic 
noise arising from low molecule numbers, as expected. 


Importantly, however, it demonstrates that when cou¬ 
pled with explicit mRNA and competitive degradation 
dynamics, these intrinsic effects dominate the response 
up to a much higher molecule number regime. 

Finally, we relax the third assumption by considering a 
three-species model for the SynExSlow circuit, in which 
the dynamics of ComS are accounted for explicitly (see 
Fig. A9). We find that effect ii persists, while effect i does 
not, indicating that the expansion of the viable response 
regime at high induction levels is more robust than at low 
induction levels. Since this is also the more pronounced 
effect, we focus on the high-induction regime in the next 
section, where we compare our model predictions with 
experiments. 


C. Fluorescence microscopy confirms the 
predictions of the model 

To test our model predictions, we use quantitative flu¬ 
orescence microscopy to measure the ComK expression 
levels in populations of B. subtilis cells harboring either 
the native or the SynExSlow circuits, as described in Sec. 
V (see Fig. 5A). ComK expression is induced by increas¬ 
ing the concentration of IPTG, which corresponds to the 
model parameter a^. As seen in Fig. 5B, in both the 
native and the SynEx strain, as the IPTG concentration 
increases, the fluorescence distribution across the popu¬ 
lation changes shape: first it is centered at low values, 
then it is split between low and high values, and finally 
it is centered at high values. This change is qualitatively 
reminiscent of the change seen in the stochastic model in 
Fig. 3B. Moreover, Fig. 5C also shows that the transi¬ 
tion to a distribution centered at high values occurs at a 
higher IPTG concentration in cells with the native circuit 
than in cells with the SynExSlow circuit. This feature is 
also qualitatively consistent with the theoretical predic¬ 
tion shown in Fig. 3C. 

To investigate whether our experimental observations 
agree quantitatively with our theoretical predictions, as 
well as qualitatively, we fit the fluorescence distributions 
to the stochastic model, as described in Sec. V. Fig. 
5G shows that both of our central predictions in the 
high-induction regime are quantitatively confirmed by 
the data, namely (i) that noise extends the transition 
to a permanently competent state beyond the determin¬ 
istically predicted induction level, and (ii) that it does so 
to a larger extent in the native circuit than in the SynEx 
circuit. Fig. 5 therefore provides strong experimental 
support for the the notion that intrinsic noise expands 
the viable response range by delaying, as a function of 
induction level, the relaxation of cells to the competent 
state. 
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FIG. 5; Quantitative fluorescence microscopy confirms model predictions. (A) Microscopy image of B. subtilis cells 
with 1.5 /iM IPTG. GomK expression, measured by CFP fluorescence, is shown in red. (B) Fluorescence distributions over the 
imaged populations for both circuits as a function of IPTG levels. Note that, as in the model, the shift from a non-responsive 
state (low fluorescence) to a responsive state (high fluorescence) is clearly evident in the distributions. (C) Comparison of 
the data with the model in the high-induction regime. For both circuits, ak is normalized by the value of the deterministic 
transition from the oscillatory to the mono-stable regime (dashed line). Agreement between the model and data confirms both 
model predictions: that noise extends the viable response regime to higher stress levels than predicted deterministically, and 
that the effect is more pronounced in the native circuit than in the SynEx circuit. 


IV. DISCUSSION 

Phenotypic heterogeneity, in which different individ¬ 
uals express particular genes at different levels, is an 
important survival strategy in uncertain environments. 
Here we studied dynamically maintained phenotypic het¬ 
erogeneity in the competence response of B. subtilis, and 


how it is influenced by intrinsic fluctuations in molecule 
numbers. By combining theoretical modeling, stochas¬ 
tic simulations, and quantitative microscopy, we showed 
that intrinsic noise facilitates heterogeneity by expand¬ 
ing the range of stress levels over which heterogeneity is 
maintained (the viable response range). The effect man¬ 
ifests itself at both low and high stress levels, and the 








































influence of noise is dramatic: in the native competence 
circuit, noise increases the maximal stress level at which 
a heterogeneous population response occurs, by 20-fold. 

Our work advances previous work investigating the 
effects of circuit architecture on dynamic response. It 
was previously known that the native competence cir¬ 
cuit exhibited higher variability in its competence du¬ 
ration times than a synthetic analog with different ar¬ 
chitecture (SynEx). This variability was attributed to 
intrinsic molecule number fluctuations and was thought 
to provide a fitness advantage, similar to variability in 
the times to commit to cell states [18]. Yet the advan¬ 
tages of the native design over the synthetic design were 
not immediately clear. Here, we showed that while the 
SynEx circuit is more predictable in terms of competence 
duration times, its dynamic response is limited to a much 
smaller range of stress levels, which limits its functional¬ 
ity. 

In both circuits, the viable response range is expanded 
at both low and high stress levels. The mechanisms in 
these two cases are different. At low stress levels, in¬ 
trinsic noise causes repeated, period excitations, effec¬ 
tively sustaining oscillations into the excitable regime. 
At high stress levels, noise prevents the damping of oscil¬ 
lations, effectively delaying, as a function of stress level, 
the static and indefinite entry into the competent state. 
Both mechanisms rely on intrinsic fluctuations and, im¬ 
portantly, persist even at high molecule numbers in a 
non-adiabatic system. The limited response range ob¬ 
served in the deterministic solution is recovered only in 
the strict limit of fast switching and high molecule num¬ 
bers, suggesting that the effects of noise that we observe 
here are generic. 

Quantitative microscopy measurements conhrmed our 
theoretical predictions: all cells exhibited competence at 
high induction levels, no cells exhibited competence at 
low induction levels, and a bimodal population response 
was observed in the intermediate regime. Important dif¬ 
ferences between circuit architectures were also conhrmed 
experimentally, namely that the SynEx circuit begins to 
oscillate at lower stress levels (IPTG levels) than the na¬ 
tive circuit, and that the native circuit can withstand 
roughly 5-fold higher stress levels than the SynEx cir¬ 
cuit before indehnitely entering the competent state. An 
independent ht of the model predictions to the experi¬ 
mental data showed very good agreement. 

Traditionally, noise in gene expression has often been 
seen as a nuisance that needs to be controlled, especially 
in stable environments or when the reproducibility of 
downstream gene expression is crucial. Thus, much work 
has concentrated on how to ensure the reliability of gene 
expression and cell signaling in the presence of intrinsic 
noise [19-27]. However, as has been shown experimen¬ 
tally and theoretically in the context of antibiotic resis¬ 
tance [3, 28, 29], noise-induced population heterogeneity 
can be advantageous for adaptation to new conditions 
[30-33]. Functional applications of noise have also been 
identihed in a number of settings [34] ranging from dif¬ 


ferentiation decisions to sporulate [35], apoptose [36], or 
allow DNA uptake, such as discussed in this paper. 

In B. subtilis, a heterogeneous competence response is 
thought to be optimal since permanent competence curbs 
cell growth. The effect of phenotypic heterogeneity on 
the growth rate of populations has also been studied the¬ 
oretically [2, 37-42], showing that while in optimal con¬ 
ditions fluctuations decrease the overall growth rate, in 
less favorable environments, diversity of gene expression 
increases the population Htness [42]. The effect of selec¬ 
tion on such populations was also considered [41], and 
shown to influence the stability of the phenotypic states 
[ 6 ]. 

We have described an example of phenotypic hetero¬ 
geneity that is maintained by an oscillatory response, and 
we have demonstrated that intrinsic noise increases the 
range of stress levels for which oscillations occur. The 
ability of noise to facilitate oscillations has also been ob¬ 
served in the entrainment of NF-kB in fibroblast cells 
to oscillating TNF inputs [43]. There, small-molecule- 
number noise was shown to facilitate both oscillation and 
entrainment, and phenotypic variability was shown to en¬ 
large the dynamic range of inputs for which entrainment 
is possible. These results, along with our findings herein, 
suggest that strategies that exploit the coupling between 
noise and phenotypic heterogeneity allow for functional 
population responses over a large variety of conditions. 

V. MATERIALS AND METHODS 
A. Stochastic model 

Our minimal stochastic models of the native and 
SynEx circuits are based on our previous modeling work 
[16], but employ the (stochastic) master equation instead 
of a (deterministic) dynamical system in order to capture 
the effects of intrinsic noise. The master equation de¬ 
scribes the dynamics of the probability distribution over 
the numbers of the relevant molecular species inside the 
cell [44]. For both circuits the master equation reads 

dpnm / I 1 \ 

— 9n—lPn—l,ni~^ ^n+l,m(^“t” l)Pn-t-l,m 

~\~QnPn,m—l 4” 'Sn.m-t-l (^ 4“ ^)Pn,m+l 

4 - Tnm'^ 4 - Qn 4 - Snm'^)Pnm' ( 1 ) 

where Pnm is the joint probability distribution over 
molecule numbers n and m (see Fig. 2 for a diagram 
and explanation of all variables and parameters). In the 
native circuit, n is the number of ComK proteins and m 
is the number of ComS proteins. In the SynEx circuit, 
n is the number of ComK proteins and m is the number 
of MecA proteins. The dynamics are birth-death pro¬ 
cesses with mutual regulation: the production rates 9 
and q increase the numbers n and m, respectively, while 
the degradation rates r and s decrease the numbers n 
and TO, respectively, and the regulation is encoded in the 



9 


functional dependence of the rates on n and m. The reg¬ 
ulation functions follow from our previous work [16] and 
for the native circuit read 


g-n — Ctk + 


qn = C(s + 


ki 


■ 


Ps 


1 -I- (n/ks)P' 

_6 k _ 

1 -t-n/Tfe + m/Ts 
6 . 


+ Afe, 


Snm — As, 

1 -I- n/1 k + m/l s 

while for the SynEx circuit they read 

PkU^ 


9n — Clk + 


ki 


Qn — “t“ 


PmnP 


km+ nP' 
Pnm — Pm — 6l7l Afc, 

^'n m. ^ Am,. 


( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 

( 7 ) 

( 8 ) 
( 9 ) 


The meaning of the parameters is explained in Fig. 2. 
The regulatory functions introduce positive and negative 
feedbacks (see Fig. 2, top). The parameter values used 
in the model are as in [16] and are given in Appendix A. 

The model in Eqns. 1-9 makes three simplifying as¬ 
sumptions, all of which we later relax. First, as in [15, 16] 
we have assumed that mRNA dynamics and the enzy¬ 
matic degradation process are substantially faster than 
all other biochemical reactions in the circuits, and are 
thus adiabatically eliminated (see Appendix A for de¬ 
tails). This reduces each model to the two-species form 
in Eqn. 1, depicted by the cartoons in Fig. 2 (top). Sec¬ 
ond, the parameters are chosen such that typical protein 
copy numbers are small (in the tens or hundreds per cell). 
Lacking information about the absolute protein numbers 
in the experiments, we make this assumption because we 
expect any effects of intrinsic noise to be most evident 
in the low-number regime, although as we later show, 
the effects we find persist out to protein numbers in the 
tens of thousands. Third, we model in Eqns. 6-9 the 
SynEx circuit as originally constructed [16], instead of 
the “SynExSlow” circuit that we use in experiments (de¬ 
scribed later in this section). This reduces the model 
from three species to two, which is more amenable to an¬ 
alytic and numerical solution. Once again, however, we 
will see that the most important effects of noise that we 
elucidate are also present in a model of the SynExSlow 
circuit. 

We solve Eqn. 1 in steady state in one of two ways. At 
low copy numbers, we use the spectral method [45, 46], 
a hybrid analytic-numerical technique that exploits the 
eigenfunctions of the birth-death process. Derivation of 
the spectral solution of Eqn. 1 is given in Appendix 
A. The spectral method is much more efficient than 
other numerical techniques [45], but we find here that 


it becomes numerically unstable at sufficiently high copy 
numbers. Therefore, at high copy numbers, we use iter¬ 
ative inversion of the matrix acting on Pnm on the right- 
hand side of Eqn. 1 (see Appendix A). To obtain indi¬ 
vidual stochastic trajectories of the system described by 
Eqns. 1-9, we use the Gillespie algorithm [47]. 


B. Deterministic model 

The deterministic analog of Eqn. 1 is obtained by per¬ 
forming an expansion in the limit of large molecule num¬ 
bers [44]. To first order one obtains 

gn - rnmfi, ( 10 ) 

Qn - Snmfh, (11) 

where h and fh are ensemble averages. Eqns. 10 and 
11 form a coupled dynamical system whose properties 
we obtain by linear stability analysis. As shown in Fig. 
Al, both circuits exhibit excitable, oscillatory, and mono¬ 
stable regimes, depending on the value of the control pa¬ 
rameter Ofc. The transition from excitable to oscillatory 
is marked by the annihilation of a stable and an unstable 
fixed point, leaving only one unstable fixed point. The 
transition from oscillatory to mono-stable is marked by 
this unstable fixed point becoming stable. These transi¬ 
tions provide the dashed lines in Fig. 3C. 


dn 

dt 

dm 

dt 


C. Genetic circuit construction 

For the native competence circuit, we used a variant 
from our previous study [13]. For the synthetic compe¬ 
tence circuit, we reconfigured the original “SynExSlow” 
circuit created in [16], in order to introduce a tunable 
proxy for stress level. This required replacing the tun¬ 
able Phyperspank ~ comS with an internally controlled 
promoter for the ribosomal gene rpsD, as well as adding 
in Phyperspank — comK. The result was a strain that 
is resistant to four antibiotics and has comS expressed 
from a ribosomal promoter, providing for a basal level of 
expression (see Appendix A for chromosomal alterations 
and antibiotic resistance). In both strains, ComK expres¬ 
sion is induced by increasing the amount of IPTG in the 
environment. Since stress signals are usually integrated 
at the ComK promoter, IPTG therefore acts as a proxy 
for stress and triggers competence. This allows us to sim¬ 
ulate stress directly in a controlled manner, rather than 
using physiological stresses that may themselves induce 
external variation in the responses. 


D. Time-lapse microscopy 

Cells of Bacillus suhtilis were prepared by streaking 
from glycerol stocks onto LB agar plates containing the 
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appropriate antibiotic for maintenance and incubated at 
37 °C overnight.. Single colonies were then selected from 
the plates and grown in LB broth for three to four hours 
at 37 °C until an OD of 1.6 to 1.8 was reached. While 
culturing the cells, argarose pads were made by pour¬ 
ing 6 mL of 0.8% w/v low-melting point agarose in re¬ 
suspension medium onto a glass coverslip. A second glass 
coverslip was then placed on top of the medium, and the 
medium was left to congeal while the culture was grown. 
Once the culture was ready, cells were spun down and 
resuspended in the resuspension medium twice to wash 
way the LB. To deposit cells, the top glass coverslip was 
removed, and then 2 /iL of cells were dropped on 37 °C 
low melting point agarose pads. The pads were then cut 
into squares with a 5mm edge, each containing a sin¬ 
gle drop of cells. After drying for one additional hour, 
the pads were flipped over and placed on a glass-bottom 
dish. The dish was then sealed with parafilm. Images of 
the cells were then obtained at lOOX magnification on an 
Olympus 1X81 system using the ImagePro software from 
MediaCybernetics along with customized macros. 

IPTG stock solutions were dissolved in ethanol to a 
concentration of 100 mM. Working (lOOOX) stocks were 
diluted with Milli-Q water to 30 mM, 10 mM, 3 mM, 1.5 
mM, 0.75 mM, respectively, by serial dilution and then 
added to the appropriate media at a ratio of 1:999 to 
achieve the final concentrations indicated. 


E. Plasmid and strain constrnction 

Template plasmids with homologous recombination 
arms for the Bacillus subtilis chromosomal loci were mod¬ 
ified through restriction enzyme digest and ligation of 
DNA inserts (see Appendix A for loci). The inserts were 
created by polymerase chain reactions using primers from 
Integrated DNA Technologies while using genomic DNA 
or other plasmids as templates. 

The PY79 strain of Bacillus subtilis was modified 
through homologous recombination using a One-Step 
Transformation protocol by inducing competence. 50 
ng of plasmid DNA was replicated in TOPIO E. coli 
cells (Invitrogen, Life Sciences, Inc) and purified using 
a MiniPrep spin column (Sigma-Aldrich). The DNA was 
then mixed with culture growing in minimal salts for 
thirty minutes and then subsequently were rescued using 
2xYT rich medium. Positive colonies were then selected 
on LB agar plates containing selective concentrations of 
antibiotics. 


F. Image analysis 

Fluorescence histograms were obtained from mi¬ 
croscopy images using a pixel-based analysis. A mask 
was created on each image to identify the areas that the 
cells occupy (see Fig. A2). A histogram of fluorescence 


intensity values was then generated for pixels within that 
area. 


G. Culture media 

Sterlini-Mandelstram Resuspension Medium was used 
during time-lapse microscopy and followed the protocol 
as in references [48, 49]. The actual protocol used con¬ 
sists of making two salt solutions: A and B. Solution A 
consists of 0.089 g of FeCl3-6H20, 0.830 g of MgCl2-6H20 
and 1.979 g MnCl 2 -4112 0 in 100 mL of filtered water. So¬ 
lution A is filter sterilized (not autoclaved) and stored at 
4 °C. Solution B consists of 53.5 g NH 4 CI, 10.6 g Na 2 S 04 , 
6.8 g KH 2 PO 4 , and 9.7 g NH 4 NO 3 . Solution B is then 
also filter sterilized and stored at 4 °C. Sporulation salts 
are made by combining adding 1 mL of Solution A and 
10 mL of Solution B to filtered water for a total 1 L. This 
solution is then autoclaved. The final Resuspension me¬ 
dia is created by combining 93 mL of sporulation salts, 2 
mL of 10% v/v L-glutamate, 1 mL of O.IM CaCl 2 , and 
4 mL of IM MgS 04 on the day of the experiment. 

One-step transformation media consists of 6.25 g of 
K 2 HPO 4 • 3 H 2 O, 1.5 g of KH 2 PO 4 , 0.25 g of trisodium 
citrate, 50 mg of MgS 04 • 7 H 2 O, 0.5 g of Na 2 S 04 at pH 
7.0, 125 of 100 mM FeCl 3 , 5 /rL of 100 mM MnS 04 , 
1 g of glucose, and 0.5 g of glutamate added into filtered 
water for a total of 250 mL. The media is filter sterilized 
using 0.2 micron Millipore filters. 

2xYT recovery medium consists of 16.0 g of Tryptone, 
10.0 g of Yeast Extract, 5.0 g of NaCl added to filtered 
water to a total volume of IL. The media is then filter 
sterilized using 0.2 micron Millipore filters. 


H. Distributiou analysis aud compariug 
experimeuts with modeliug 

For the ComK distributions in the model, = 
J2mPnm, we determine the fraction / of in the respon¬ 
sive, high comK protein concentration state using two in¬ 
dependent methods. First, we use a generalized method 
of separating the distribution’s two modes: since the dis¬ 
tribution is often not completely bimodal (see Fig. 3B, 
middle column), we find the average n* = {rii + n 2 )f 2 
of the two inflection points surrounding the putative lo¬ 
cal minimum between the two modes, and define / = 
Pn- III Ih6 case of a bimodal this method in¬ 
deed well approximates the location of the actual local 
minimum. Second, we fit to a mixture of two Pois¬ 
son distributions, using the Kullback-Leibler divergence 
as the cost function. There are three fitting parameters, 
the two Poisson parameters and the relative weighting be¬ 
tween them, and the weighting provides /. We see in Fig. 
A3 that the two methods give similar results for the de¬ 
pendence of / on the control parameter a^, demonstrat¬ 
ing that our determination of / is robust to the method 
used. 
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To compare the model predictions to the experimental 
data, we analyze the fluorescence distributions (Fig. 5B) 
in the same way as the model distributions. Specifically, 
we calculate the fraction of the comK population in the 
responsive state / for each experimental distribution us¬ 
ing the two-Poisson method above. Because the mapping 
between IPTG concentration and the model parameter 
describing the level of external stress Uk is unknown, we 
infer the most likely value of corresponding to each ex¬ 
perimental distribution by fitting the theory to the data. 
First we use the mode of the [IPTG] = 0 distributions to 
subtract the background fluorescence from the remaining 
data. To avoid binning, we then fit the cumulative dis¬ 
tribution instead of the probability distribution (sample 
fits are shown in Fig. A4). We use a maximally con¬ 
strained least-squares fit, where all parameters are fixed 
as in Appendix A except and the unknown parameter 
X describing the conversion of pixel intensity to molecule 
number. Given a value of X, we find the values of ak for 
each distribution that minimize the sum of the squared 
error S in each case. X is then chosen by minimizing the 


sum of minimum S values over all distributions. These 
/ and ak values inferred from the data are plotted in 
Fig. 5C. The error bars on ak are obtained by finding 
the ak values where S reaches 1.25 of its minimum value 
(sample plots of S vs. ak are shown in Fig. A4). 


Acknowledgments 

M.K. was supported by NIH biophysics training grant 
T32GM008297 and acknowledges Tolga Cagatay for as¬ 
sistance in preparing Bacillus strains. L.H. was sup¬ 
ported by NSF REU grant PHY-1460899. G.M.S. was 
supported by the National Institute of General Medi¬ 
cal Sciences Grant ROl GM088428 and the National Sci¬ 
ence Foundation Grant MCB-1450867. A.M.W. was sup¬ 
ported by MGGIG grant no. 303561. This work was also 
supported by the San Diego Center for Systems Biology 
(NIH Grant P50 GM085764). 


[1] Dawn Fraser and Mads Kaern. A chance at survival: gene 
expression noise and phenotypic diversification strate¬ 
gies. Molecular microbiology, 71 (6): 1333-1340, 2009. 

[2] Edo Kussell and Stanislas Leibler. Phenotypic diversity, 
population growth, and information in fluctuating envi¬ 
ronments. Science, 309(5743):2075-2078, 2005. 

[3] Nathalie Q Balaban, Jack Merrin, Remy Chait, Lukasz 
Kowalik, and Stanislas Leibler. Bacterial persistence as a 
phenotypic switch. Science, 305(5690) :1622-1625, 2004. 

[4] Anna Kuchina, Lorena Espinar, Jordi Garcia-Ojalvo, and 
GM Suel. Reversible and noisy progression towards a 
commitment point enables adaptable and reliable cellu¬ 
lar decision-making. PLoS Comput Biol, 7(ll):el002273, 
2011 . 

[5] Anna Kuchina, Lorena Espinar, Tolga Qagatay, Ale¬ 
jandro O Balbin, Fang Zhang, Alma Alvarado, Jordi 
Garcia-Ojalvo, and Giirol M Siiel. Temporal competi¬ 
tion between differentiation programs determines cell fate 
choice. Molecular systems biology, 7(1):557, 2011. 

[6] Thierry Mora and Aleksandra M Walczak. Effect of 
phenotypic selection on stochastic gene expression. The 
Journal of Physical chemistry B, 117(42): 13194-13205, 
2013. 

[7] Aaron Novick and Milton Weiner. Enzyme induction as 
an all-or-none phenomenon. Proceedings of the National 
Academy of Sciences of the United States of America, 
43(7):553, 1957. 

[8] Paul J Ghoi, Long Gai, Kirsten Frieda, and X Sunney 
Xie. A stochastic single-molecule event triggers pheno¬ 
type switching of a bacterial cell. Science, 322 (5900) :442- 
446, 2008. 

[9] Michael A Savageau. Gomparison of classical and autoge¬ 
nous systems of regulation in inducible operons. Nature, 
252:546-549, 1974. 

[10] Mark Kittisopikul and Giirol M Siiel. Biological role of 
noise encoded in a genetic network motif. Proceedings of 
the National Academy of Sciences, 107(30): 13300-13305, 


2010 . 

[11] Rachel S Koh and Mary J Dunlop. Modeling suggests 
that gene circuit architecture controls phenotypic vari¬ 
ability in a bacterial persistence network. BMC systems 
biology, 6(1):47, 2012. 

[12] Denis Michel. Kinetic approaches to lactose operon in¬ 
duction and bimodality. Journal of theoretical biology, 
325:62-75, 2013. 

[13] Giirol M Siiel, Rajan P Kulkarni, Jonathan Dworkin, 
Jordi Garcia-Ojalvo, and Michael B Elowitz. Tunability 
and noise dependence in differentiation dynamics. Sci¬ 
ence, 315(5819):1716-1719, 2007. 

[14] B-J Haijema, J Hahn, J Haynes, and D Dubnau. A 
comga-dependent checkpoint limits growth during the es¬ 
cape from competence. Molecular microbiology, 40(1):52- 
64, 2001. 

[15] Giirol M Siiel, Jordi Garcia-Ojalvo, Louisa M Liberman, 
and Michael B Elowitz. An excitable gene regulatory 
circuit induces transient cellular differentiation. Nature, 
440(7083):545-550, 2006. 

[16] Tolga Qagatay, Marc Turcotte, Michael B Elowitz, 
Jordi Garcia-Ojalvo, and Giirol M Siiel. Architecture- 
dependent noise discriminates functionally analogous dif¬ 
ferentiation circuits. Cell, 139(3):512-522, 2009. 

[17] Andrew J Black and Alan J McKane. Stochastic formu¬ 
lation of ecological models and their applications. Trends 
in ecology & evolution, 27(6):337-345, 2012. 

[18] Iftach Nachman, Aviv Regev, and Sharad Ramanathan. 
Dissecting timing variability in yeast meiosis. Cell, 
131(3):544-556, 2007. 

[19] Andrew Mugler, Filipe Tostevin, and Pieter Rein ten 
Wolde. Spatial partitioning improves the reliability 
of biochemical signaling. Proceedings of the National 
Academy of Sciences, 110(15):5927-5932, 2013. 

[20] Christopher C Govern and Pieter Rein ten Wolde. 
Optimal resource allocation in cellular sensing sys¬ 
tems. Proceedings of the National Academy of Sciences, 



12 


111(49):17486-17491, 2014. 

[21] Christopher C Govern and Pieter Rein ten Wolde. En¬ 
ergy dissipation and noise correlations in biochemical 
sensing. Physical review letters, 113(25):258102, 2014. 

[22] Filipe Tostevin, Pieter Reinten Wolde, and Martin 
Howard. Fundamental limits to position determination 
by concentration gradients. PLoS Computational Biol¬ 
ogy, 3(4), 2007. 

[23] Timothy E Saunders and Martin Howard. Morphogen 
profiles can be optimized to buffer against noise. Physical 
Review E, 80(4):41902, 2009. 

[24] Gasper Tkacik, Aleksandra M Walczak, and William 
Bialek. Optimizing information flow in small genetic net¬ 
works. Physical Review E, 80(3):031920, 2009. 

[25] Aleksandra M Walczak, Gasper Tkacik, and William 
Bialek. Optimizing information flow in small genetic net¬ 
works. ii. feed-forward interactions. Physical Review E, 
81(4):041905, 2010. 

[26] Gasper Tkacik, Aleksandra M Walczak, and William 
Bialek. Optimizing information flow in small genetic 
networks, iii. a self-interacting gene. Physical Review E, 
85(4):041903, 2012. 

[27] Gasper Tkacik and Aleksandra M Walczak. Information 
transmission in genetic regulatory networks: a review. 
Journal of Physics: Condensed Matter, 23(15):153102, 
2011 . 

[28] Eitan Rotem, Adiel Loinger, Irine Ronin, Irit Levin- 
Reisman, Ghana Gabay, Noam Shoresh, Ofer Biham, and 
Nathalie Q Balaban. Regulation of phenotypic variabil¬ 
ity by a threshold-based mechanism underlies bacterial 
persistence. Proceedings of the National Academy of Sci¬ 
ences, 107(28):12541-12546, 2010. 

[29] Orit Gefen and Nathalie Q Balaban. The importance 
of being persistent: heterogeneity of bacterial popula¬ 
tions under antibiotic stress. EEMS microbiology reviews, 
33(4):704-717, 2009. 

[30] Katsuhiko Sato, Yoichiro Ito, Tetsuya Yomo, and Kuni- 
hiko Kaneko. On the relation between fluctuation and 
response in biological systems. Proceedings of the Na¬ 
tional Academy of Sciences, 100(24):14086-14090, 2003. 

[31] Yoichiro Ito, Hitoshi Toyota, Kunihiko Kaneko, and Tet¬ 
suya Yomo. How selection affects phenotypic fluctuation. 
Mol Syst Biol, 5:264, 2009. 

[32] Akiko Kashiwagi, Itaru Urabe, Kunihiko Kaneko, and 
Tetsuya Yomo. Adaptive response of a gene network to 
environmental changes by fitness-induced attractor selec¬ 
tion. PloS one, l(l):e49, 2006. 

[33] Y Shimizu, S Tsuru, Y Ito, BW Ying, T Yomo, and Grze- 
gorz Kudla. Stochastic switching induced adaptation in 
a starved escherichia coli population. PLoS ONE, 6(9), 
2011 . 

[34] Avigdor Eldar and Michael B Elowitz. Functional roles 
for noise in genetic circuits. Nature, 467(7312):167-173, 
2010 . 


[35] Jan-Willem Veening, Eric J Stewart, Thomas W Bern- 
gruber, Frangois Taddei, Oscar P Kuipers, and Leen- 
dert W Hamoen. Bet-hedging and epigenetic inheritance 
in bacterial cell development. Proceedings of the National 
Academy of Sciences, 105(ll):4393-4398, 2008. 

[36] Sabrina L Spencer, Suzanne Gaudet, John G Albeck, 
John M Burke, and Peter K Sorger. Non-genetic ori¬ 
gins of cell-to-cell variability in trail-induced apoptosis. 
Nature, 459(7245):428-432, 2009. 

[37] Mukund Thattai and Alexander Van Oudenaarden. 
Stochastic gene expression in fluctuating environments. 
Cenetics, 167(l):523-530, 2004. 

[38] Edo Kussell, Roy Kishony, Nathalie Q Balaban, and 
Stanislas Leibler. Bacterial persistence a model of sur¬ 
vival in changing environments. Genetics, 169(4): 1807- 
1814, 2005. 

[39] Stanislas Leibler and Edo Kussell. Individual histories 
and selection in heterogeneous populations. Proceed¬ 
ings of the National Academy of Sciences, 107(29): 13183- 
13188, 2010. 

[40] Olivier Rivoire and Stanislas Leibler. The value of infor¬ 
mation for populations in varying environments. Journal 
of Statistical Physics, 142(6):1124-1166, 2011. 

[41] Katsuhiko Sato and Kunihiko Kaneko. On the distribu¬ 
tion of state values of reproducing cells. Physical biology, 
3(1):74, 2006. 

[42] Sorin Tanase-Nicola and Pieter Rein Ten Wolde. Regu¬ 
latory control and the costs and benefits of biochemical 
noise. PLoS Comput Biol, 4(8):el000125-el000125, 2008. 

[43] Ryan A Kellogg and Sava§ Tay. Noise facilitates 
transcriptional control under dynamic inputs. Cell, 
160(3):381-392, 2015. 

[44] Nicolaas Godfried Van Kampen. Stochastic processes in 
physics and chemistry, volume 1. Elsevier, 1992. 

[45] Aleksandra M Walczak, Andrew Mugler, and Chris H 
Wiggins. A stochastic spectral analysis of transcrip¬ 
tional regulatory cascades. Proceedings of the National 
Academy of Sciences, 106(16):6529-6534, 2009. 

[46] Andrew Mugler, Aleksandra M Walczak, and Chris H 
Wiggins. Spectral solutions to stochastic models of gene 
expression with bursts and regulation. Physical Review 
E, 80(4):041921, 2009. 

[47] Daniel T Gillespie. Exact stochastic simulation of cou¬ 
pled chemical reactions. The journal of physical chem¬ 
istry, 81(25):2340-2361, 1977. 

[48] J. M. Sterlini and J. Mandelstam. Commitment to sporu- 
lation in bacillus subtilis and its relationship to develop¬ 
ment of actinomycin resistance. Biochem J, 113(l):29-37, 
Jun 1969. 

[49] C.R. Harwood and S.M. Cutting. Molecular Biological 
Methods for Bacillus. John Wiley & Sons Ltd, West Sus¬ 
sex, England, 1990. 


Appendix A: Supporting Information 
1. Parameter values 


The parameter values used in the model in the main text are taken from [13, 15, 16], where they have been 
optimized to agree with a number of experimental observations, including the existence of competence events, the 
duration of these events, and the robustness of events to a wide range of stress levels. They are first given below in 
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dimensionless form, which is sufficient to establish their deterministic dynamical properties. To account for intrinsic 
noise and to compare with experiments, several physical quantities are then specified to establish molecule numbers 
and timescales. The conversion to the remaining model parameters, and their resulting values, are then given below. 


Native circuit: 


Dimensionless parameters 

Physical quantities 

Remaining parameters 

Og = 0 

Molecule numbers: 

ak (varied) 

bk = 0.3 

Pfc = 100 

ttg = OgTgd = 0/s 

=3 

rg = i 

/3k = bkTkSk = 0.03/s 

fco = 0.2 


/3g = 6grg4 = 0.003/s 

ki = 1/30 

Timescales: 

kk = koVk = 20 

Afe = 0.1 

4 = 0.001/s 

ks = kiTk = 3.3 

Ag =0.1 

4 = 0.001/s 

Afc = AkSk = 10-4/8 

h = 2 


Ag = Ag4 = 10-4/s 

p = 5 




SynEx circuit: 


Dimensionless parameters 

Physical quantities 

Remaining parameters 

— 0-3 

Molecule numbers: 

ak (varied) 

bk = 15 

kk = 10 

~ 1-5 X 10 /s 

bm = 10 


I3k = bkkkXk = 0.015/s 

p = 1 


i^m — — 0.005/s 

h = 2 

Timescales: 

d = pXk/km = 2 X 10-^/s 

p = 2 

Afe = 10-4/s 

Xm = 10-4/s 



Note that in [13, 15, 16], the higher values Tj, = 25000 and Tg = 20 (native), and kk = 5000 and km = 2500 (SynEx) 
are used to establish molecule number. We use the lower values here to elucidate the effects of intrinsic noise. This 
assumption is relaxed in the next section, where we return to the high-number regime. In the main text, we discuss 
how the effects of noise are robust to this choice. 


2. Relaxing the model assumptions 


To relax the first two simplifying assumptions made in the main text, we consider a stochastic model that has 
high molecular copy numbers, and that accounts explicitly for the mRNA dynamics and the dynamics of competitive 
degradation. The model follows from our earlier work [16], and consists of a set of coupled chemical reactions for each 
circuit. Due to the complexity of the model, we do not solve the master equation explicitly, but rather we simulate 
the dynamics using the Gillespie algorithm [47]. 


For the native circuit, the reactions and rates are: 
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Reaction 

Rate 

Additional parameters 

PcomK ^ PcomK + ^RNAcon^K 

ki (varied) 

kg = 0.19/s 

PComK —> PComK + HlRNAcomK 

fc2n^/(fc(( + n^) 

kz = 0.0015/s 

rnRNAcomK ^ mRNAconiK + ComK 

ks = 0.2/s 

fell = 2 X 10-Vs 

P&°o"^‘s ^ PcZs + mRNAcomS 

ki = 0/s 

ki3 = 4.5 X 10-Vs 

PComS —> PcomS + mRNAcornS 

fcs/]! + [n/ksY] 

kk = 5000 

mRNAcomS -> mRNAcornS + ComS 

ke = 0.2/s 

ks = 833 

mRNAcornK 0 

fcr = 0.005/s 

h = 2 

ComK —> 0 

ks = 10“^/s 

p = 5 

mRNAcornS 0 

kg = 0.005/s 

H = 1.66 pm^ 

ComS —)■ 0 

fcio = 10"'‘/s 

Mt = 500 

MecA + ComK —>■ MecAx 

kii/9l 


MccAk —>■ MecA + ComK 

k-ii = 5 X 10-Vs 


MccAk —>■ MecA 

fci 2 = 0.05/s 


MecA + ComS —)■ MecAs 

kis/fl 


MecAs —)■ MecA + ComS 

fc_i3 = 5 X 10-Vs 


MecAs —)■ MecA 

fci4 = 4 X 10-^/s 



Here Pgene denotes the promoter of the corresponding gene (constitutive or regulated), and MecAx and MecAs 
represent the complex of MecA bound to ComK and ComS, respectively. As in the main text, n is the number of 
ComK molecules per cell. H represents the cell volume, and Mr gives the total number of MecA molecules. Upon 
adiabatically eliminating the faster mRNA dynamics and the MecA dynamics, one obtains the model in the previous 
section, except with the larger molecule numbers = 25000 and Fg = 20 [13]. The control parameter here, ki, is 
related to the control parameter in the reduced model, ak, via fci = k^ak/ks [13]. 

For the SynEx circuit, the reactions and rates are: 


Reaction 

Rate 

Additional parameters 

PcoVk ^ PcomK + mRNAcomK 

ki (varied) 

kg = 0.19/s 

PComK —> PcomK + mPNAcornK 

kgrY/{k^ + rY) 

ki = 0.0625/s 

PSlVcl ^ PS + mRNAMecA 

kg = 0.0019/s 

fcr = 4 X 10-Vs 

PmbcA PmocA + hiRNAmocA 

ki-nPUkPm + nP) 

kk = 5000 

mRNAcornK -> mRNAcornK + ComK 

kg = 0.2/s 

km = 2500 

mRNAMecA —>■ mRNAMecA + MecA 

kg = 0.2/s 

h = 2 

MecA + ComK —> MecA 

ky/V, 

p = 2 

mRNAcornK 0 

ks = 0.005/s 

n = 1.66 iiuY 

mRNAMecA 0 

kg = 0.005/s 


MecA —)• 0 

kio = lO-'^/s 


ComK —> 0 

fell = 10-^/s 



Once again, Fgene denotes the promoter of the corresponding gene, n is the number of ComK molecules per cell, and 
H represents the cell volume. Upon adiabatically eliminating the faster mRNA dynamics, one obtains the model in 
the previous section, except with the larger molecule numbers kk = 5000 and km = 2500 [16]. The control parameter 
here, ki, is related to the control parameter in the reduced model, Ofc, via ki = k^au/k^. 

To relax the third simplifying assumption made in the main text, we consider a stochastic model of the SynExSlow 
circuit. The model follows from our earlier work [16], and is similar to the stochastic model in the main text, except 
that there are three species with dynamic molecule numbers: ComK (n), MecA (m), and ComS {£). The production 
rate functions are, respectively. 


9n — 


Pkn^ 

k^+n’^' 


Qn — 


PmnP 
km+nP' 


, Psn’^ 


(Al) 
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and the degradation rate functions are, respectively, 
_ 5km ^ ^ 

Tnmt — /-p I /) /-n 

1 + n/rfc + £/rs 


The parameter values are [16]: 


Parameters 

O^k 

(varied) 

13s 

= 0.5/s 


= 0.075/s 

Tfc 

= 25000 

CX.S 

= 0.5/s 

r« 

= 20 

5k - 

= 6,=2x lO-s/s 

kk 

= 5000 

Afe 

= A™ = A, = 10-Vs 

km 

= 2500 

/3k 

= 7.5/s 

ks 

= 500 

/dm 

= 2.5/s 

h = 

--p = 2 


'^nml 


6 sm 

1 + n/Tk + e/T, 


+ Ag. 


(A2) 


Note that these values correspond to the high-molecule-number regime of the native and SynEx models of the 
main text, and thus this model of the SynExSlow circuit also relaxes the first simplifying assumption of low molecule 
number. 

Finally, we note that at these parameter values, the deterministic analog of this model predicts that as a function 
of the control parameter a^, the excitable regime transitions directly into a damped oscillatory (mono-stable) regime, 
where two of the three eigenvalues of the Jacobian matrix are complex, and all have negative real part. That is, 
there is no standard oscillatory regime. Therefore, we define a heuristic boundary a),' = 0.15/s after which the 
damping is clearly evident within the first 24 hours of the deterministic dynamics (see Fig. S7). The fact that the 
stochastic dynamics exhibit sustained oscillations in the absence of a deterministically oscillatory regime, even beyond 
this heuristic boundary (Fig. S7), indicates that the ability of noise to induce oscillations in the SynExSlow circuit is 
especially strong. 


3. Spectral solution to the master equation 

We write the master equation (Eqn. 1) as 


dp„ 


dt 


(Cn [9nj '^nm] [Qu} ^nm]^ Pnm- 


(A3) 


Here —C is the linear birth-death operator, whose action on the probability distribution p is described by 


-jOn[gn,rn]Pn = g-n-lPn-l + ?’n+l(n + ^)Pn+l “ (ffn + r„n)p„, (A4) 

where in general on an operator (here, C) we use a subscript to denote the sector (n or m) on which it acts. Explicitly, 
then, the master equation reads 

gn—lPn—l,m “f ^n-t-l,m (?1 H- l)j5n+l,m i.9n nm 

~^gnPn,m—l “t” ^n,m+l{m J- l')Pn^m+l {Qu T ^nmm^Pnm- (-^-5) 

We will now derive the spectral decomposition of the master equation by introducing the generating function. For 
intuition, we will first introduce the generating function in the context of the one-dimensional system described by 
Eqn. A4, then extend our results to the full master equation. 

The generating function is an expansion in a complete set of states, indexed by molecule number, for which the 
probabilities provide the expansion coefficients. Denoting the states abstractly as |n), the generating function is 
defined 



dt 


|G) = ^p„|n). 


n 


(A6) 
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where the sum runs from 0 to oo (as do all sums hereafter unless otherwise specified). Summing Eqn. A4 against \n) 
yields 

-t\G) = ^g„-iP„-i|n) + ^r„+i(n + l)p„+i|n) - ^ 5 „p„|n) - ^r„np„|n) (A7) 

n n n n 

= '^9nPn\n+l)+'^rnnpn\n-l)-'^gnPn\n)-'^rnUPnlri), (A8) 


where the second step shifts the sum without penalty (i) in the first term, since p_i =0, and (ii) in the second term, 
since it vanishes for n = 0. Eqn. A8 benefits from the introduction of two additional sets of operators: (i) operators 
corresponding to the evaluation of the production and degradation functions at particular values of n. 


g\n) = 5 „|n), (A9) 

r|n) = r„|n), (AlO) 

and (ii) raising and lowering operators that correspond to the shifts in n caused by birth and death, 

d+|n) = |n + l), (nla"*" = (n — 1|, (All) 

a“|n) = n|n — 1), (n|d“ = n(n + 1|. (A12) 

In Eqns. A11-A12, for completeness, we have also presented the operators’ actions to the left, which derive from the 
orthonormality condition (n|n') = 5nn'- In terms of the above operators, Eqn. A8 becomes 

-C\G) = ^p„d+g|n) + ^p„d“r|n) - ^p„5|n) - ^p„d+d“r|n), (A13) 

n n n n 

= {a''~g + a~f - g - a''~a~r)'^pn\n), (A14) 

n 

= - {a+-l) {a~f - g)\G). (A15) 

Eqn. A15 reveals the form of the birth-death operator in generating function space. 


We now generalize Eqn. A15 to the two-dimensional master equation of the model. Defining the two-dimensional 
generating function. 


|G) = y^Pnm|n,m), 

nm 


(A16) 


the master equation becomes 

|G) = ~ [(^n ~ l) (p'n^nm ~ Pn) + ~ l) ~ 9n)] |l^)) (A17) 

where dot denotes the time derivative, and as before the subscripts on operators denote the sectors on which they 
act. 

The key insight of the spectral method is that a master equation such as Eqn. A17 can be simplified significantly 
by expansion in a wisely chosen set of eigenfunctions. Since Eqn. A17 describes two coupled birth-death processes, 
we choose to expand in the eigenfunctions of two uncoupled birth-death processes - that is, processes with constant 
production rates g and q and degradation rates f and s, respectively. Again for intuition we begin in one dimension, 
for which the operator describing a constant-rate birth-death process is given by Eqn. A15: 

-C\G) = - (a+ - 1) (d-r-g) |G) = -fb+b-\G). (A18) 

In the second step we have defined the shifted raising and lowering operators 

b+ = d+ - 1, (A19) 

b~ = d~ — g/f. (A20) 

The operator b'^b~ is a number operator whose eigenvalues are integers, which we call j; thus the eigenvalue relation 
for the constant-rate birth-death operator is 


■^\j) = -rj\j)- 


(A21) 
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Accordingly, &+ and b raise and lower the eigenstates as a+ and a do the |n) states (Eqn. A11-A12): 

= |j + 1). = 




{j\b =j(j + l|- 


(A22) 

(A23) 


As we will see, the spectral method exploits the expansion of the generating function in these eigenstates \j). 

Returning to two dimensions, it is clear that we would like to write the coupled master equation in terms of the 
uncoupled rates, which we do by adding them to and subtracting them from Eqn. A17: 

|G) = - [(a+ - 1) {a~[fnm + r - f] - - 5 + 5 ) + (a+ - 1) + s - s] - qn - q + q)] jG) (A24) 

(A25) 

(A26) 


1 ) ^ 9 + (^m l) 9 + A„^ |G) 


rb^b — b^ 


+ f) 


+1) 


^nm 


\G). 


7nm + ~ 

In the second step we define operators which capture the deviations between the constants and the coupled rates, 

Tn = 9 9n'j = Q 9m (A27) 

7nm = T Tnm? ^nni = S Snrm (A28) 

and in the third step we write the raising and lowering operators in terms of their shifted counterparts: 

b+ = a+- 1, b+=a+- 1, (A29) 

bn=a~-9lr, Kn = Kn-9ls. (A30) 

We now expand the generating function in the eigenstates of b^b~ and 6+5“, with eigenvalues j and k, respectively, 

(A31) 


\G) = J2Gjk\j,k), 

jk 


which makes Eqn. A26 


^ ^ Gjk ll; k) — ^ ^ Gjk b^b^'jYim _5yjynm A 5^Eyj + s/c b^b^Xjun -b^Xnm b^A^nj \j:k). (A32) 

jk jk 


Projecting from the left with the state (j, k\, 

G,k = -{fj + sk)Gjk - E Gyfe(j|5+f„|j') - E 0’ k\b+K\f, k') 


j'k' 


+ X] Gj'k'{j, fc|5+5„7„m|/, C) + I ^ Gyfc'(j, k\bl^nm\j', k') 


J'k' 


j'k' 


+ ^Gyfe,(j,fc|5+5-A„m|/,C) + |^Gyfc.(j,fc|5+A„™|/,C), 

j'k' j'k' 


(A33) 


and noting the actions of 5+ and 5 to the left (Eqns. A29-A30) yields the dynamics of the expansion coefficients, 

Gjk = -{fj + sk)Gjk - ^ [Tj^ijiGj'k + ^jj'Gj'^k-i] + ^ 

j' j'k' 

where the deviations, rotated into eigenspace, have become matrices and tensors: 


^jj' = (j|r„|/) = '^{j\n){g - gn){n\jj, 

Kjj' = (j|A„|/) = ^(i|n)(g- g„)(n|/). 

(A35) 

ifj'' = {hk\lnm\3' ,kj 

X)f, = (j,A:|A„m|/,C) 


nm 

= y^jj\n){k\m){s - Srim){n\jj{m\kj. 
nm 

(A36) 




-I- 


l.\kk' I 9 ^k-l.k' 


Gj'k' 


(A34) 
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With Eqn. A34, we have turned the original master equation (Eqn. A17) into a linear algebraic equation, involving 
only matrix and tensor multiplication. In steady state Eqn. A34 can be solved perturbatively by treating the first 
term on the right-hand side as larger than the others and iterating to convergence. 

Algorithmically, then, the procedure prescribed by the spectral method is the following. First we compute the 
deviation matrices and tensors from the given rate functions (/„, Qn, and Snm, and our chosen ‘gauges’ g, r, q, 
and s. Then we solve for the steady state of Eqn. A34 by iteration to obtain the expansion coefficients Gjk- Finally, 
we compute the probability distribution by taking the inverse transform: 

Pnm — E Gjk{n\j){m\k). (A37) 

jk 


The probability distribution provides the complete stochastic description of the steady-state process. Note that both 
the deviation matrices/tensors and the probability distribution are computed via multiplication against the eigenmodes 
(n|j) and {m\k) or their conjugates (jjn) and {k\m); these are efficiently precomputed via recursive update rules. 

For the SynEx circuit, the dynamics of the spectral expansion coefficients are simplified because we have —>■ 
and Snm s: 


Gjk — -{rj -I- sk)Gjk - ^ [Tj-ijiGj'k + ^jj'Gj'^k-i] + 

j' k 


E[. 


Jlkk'Gjk' + j-^kk'Gj-i^k' 


(A38) 


where 


- 9n){n\j'), 
n 

Ikk' = ^(fc|m)(r - rm){m\k'). 


- <ln){n\j'), (A39) 

n 

(A40) 


We see that the sparser coupling of the degradation terms (compared to the native circuit) has left us with a linear 
algebraic equation involving only matrices, not tensors. The solution follows the steps outlined above for the native 
circuit. 

The spectral method is more efficient than other solution techniques by many orders of magnitude [45]. However, 
we find that here it becomes numerically unstable at sufficiently high numbers. Therefore, at high copy numbers, 
we solve for the steady state of the original master equation (Eqn. A3) by iteration. This is less efficient but more 
numerically stable. Specifically, we tile the columns of Pnm into one vector Pi, and write the master equation as 


dP^ 

dt 




{—DiSw + Nii') Pi', 


(A41) 


where Mu' is the reshaped operator in Eqn. A3, and —Di and Nn' are its diagonal and non-diagonal components, 
respectively. In steady state, we have 


i' 


(A42) 


which we solve iteratively until convergence. 


4. Chromosomal alterations and antibiotic resistance 

The native strain was a variant from a previous study [13] and had the following chromosomal alterations and 
antibiotic resistance: 


Locus 

Construct 

Antibiotic Resistance 

AmyE 

Phyperspank COmK 

Spectinomycin 

SacA 

PcomG cfp, PcomS Y^P 

Chloramphenicol 


The SynExSlow strain was modified from [16] as described in the main text and had the following chromosomal 
alterations and antibiotic resistance: 
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Locus 

Construct 

Antibiotic Resistance 

AmyE 

Phyper 

Spectinomycin 

Sac A 

PcomGKboxl'kklGcA. ^ 

Chloramphenicol 

Asr/A, comS 

PcomG^fP 

Neomycin / Kanamycin 

GltA 

PrpsDCOTTlS^ PcomG^^^^^ 

Phleomycin 


5. Supplementary figures 




ComK induction rate, ak (hr 


FIG. Al: Linear stability analysis of deterministic system predicts three dynamic regimes, for both the native 
and SynEx circuits. Fixed points h* satisfying the steady state of Eqns. 10 and 11 for (A) the native circuit and (B) the 
SynEx circuit. Each fixed point is stable if the real parts of the eigenvalues of the Jacobian matrix evaluated at that point are 
negative, and unstable otherwise. The Jacobian matrix is Jij = dFi/dxj, where x\ = n, X 2 = rfi, F\ = dn/dt, and F 2 = dfh/dt. 
For both circuits, there are three dynamic regimes. The excitable regime (low Uk) is has three fixed points, one of which is 
stable. The oscillatory regime (intermediate ak) has one unstable fixed point. The mono-stable regime (high ak) has one stable 
fixed point. In the mono-stable regime, near the oscillatory regime, the eigenvalues are complex, indicating damped oscillations. 
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1. Filter connected areas based on region properties 




3. Remove coloring to get binary mask 


2 . 





4. Apply binary mask to fluorescence image to 
generate pixel based histogram 



Contiguous areas are colored the same 


FIG. A2: Image analysis procedure. Fluorescence histograms are generated via isolating cells from the background by 
identifying connected and contiguous areas, applying a binary mask, and binning the resulting pixel intensities. Only pixels 
that fall within cells, not the background, are included. 
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ComK induction rate, ak (normalized) 


FIG. A3: Determination of high-response fraction / is robust to calculation method. Two independent methods 
are used to determine /: finding the inflection points, and fitting to a mixture of two Poisson distributions (see Sec. V). For 
both the native and SynEx circuit, we see that the two methods give results that correspond very closely to each other. The 
two-Poisson method gives smoother results, but does not capture the transition to / = 1, since a roughly equal mixture of two 
Poisson distributions with similar means (/ ~ 0.5) will always provide a better fit than a single Poisson distribution (/ = 1). 
Therefore, in Fig. 3C and Fig. 5C, we use the two-Poisson method for / < 1 and use the inflection-point method to determine 
the transition to / = 1. 


Native 


SynEx 







FIG. A4: Sample fits to experimental data and plots of sum-of-squares. Gumulative probability distributions of 
fluorescence data are fit to the theoretical distributions by minimizing the sum of squared errors, for the (A-D) native and 
(E-H) SynEx circuits, as described in Sec. V. 
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ki = 0.2 X k\‘ 


,(i) 


ki = 0.9 X k\' 


,(i) 





0 0.05 0.1 0 0.05 0.1 

frequency (mHz) frequency (mHz) 




0 0.05 0.1 

frequency (mHz) 



ki = 50 X Ari' 


( 2 ) 


-r 2 



100 

Time (hours) 


200 0 


100 

Time (hours) 


200 





0 0.05 0.1 0 0.05 0.1 

frequency (mHz) frequency (mHz) 


FIG. A5: Effects of noise persist when assumptions are relaxed in the model of the native circuit. Top two 

rows show ComK and ComS time series from Gillespie simulations of the relaxed model that includes mRNA and competitive 
degradation dynamics (see Sec. A 2). Far from the deterministic boundaries of the control parameter, and (indicated 
by the dashed vertical lines), the dynamics are excitable, oscillatory, and mono-stable as predicted (columns 1, 3, and 5, 
respectively). However, near the boundaries, but outside the oscillatory regime, noise causes oscillations to persist, due to 
either repeated excitations (column 2) or prevention of damping (column 4), confirming the effects seen in the reduced model 
of the main text. The persistence of oscillations is verified by computing the power spectrum P{uj) = |h(u))p from the Fourier 
transform of the GomK time series n{t). For periodic signals, the power spectrum is peaked at a non-zero frequency to (and in 
some cases its harmonics). Red line is a Gaussian fit to aid the eye. 
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h = 0.9 X Ai' 


,(i) 


Time (hours) 


xIO® 


5 

o o 
Q. ^ 


L. 
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ki = yj X fcf ^ 
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Time (hours) 


xIO" 
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FIG. A6: Effects of noise persist when assumptions are relaxed in the model of the SynEx circuit. As in Fig. 
A5 but for the SynEx circuit. Once again, oscillations persist outside the deterministic boundaries as indicated by the peaked 
power spectra. Note, however, that oscillations are damped at the value ki = 5k{ ' here, whereas in the native circuit they 
persist beyond this value (see Fig. A5). This confirms the effect seen in the main text that the prevention of damping is more 
pronounced in the native circuit than in the SynEx circuit. 
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FIG. A7: Effects of noise diminish when molecule number is raised in the adiabatic model of the native circuit. 

Gillespie simulations of the adiabatically reduced model of the native circuit (as in Fig. 4), but for high molecule numbers 
(Ffe = 25000 and Fa = 20). We see that 10% outside the deterministically oscillatory regime, the stochastic dynamics are either 
non-oscillatory (column 2) or weakly oscillatory (column 4). This is in contrast to the low-molecule-number regime (Fig. 4), 
where oscillations persist in these regions and beyond. We conclude that raising molecule number in the adiabatic model of 
the native circuit reduces the stochastic behavior to the deterministic behavior. 
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cKfc = 0.9 X 


Power spectrum 
I Gaussian fit 
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FIG. A8: Effects of noise diminish when molecule number is raised in the adiabatic model of the SynEx circuit. 

As in Fig. A7, but for the adiabatically reduced model of the SynEx circuit at high molecule numbers {kk = 5000 and 
km = 2500). Comparing to Fig. 4, we similarly conclude that raising molecule number in the adiabatic model of the SynEx 
circuit reduces the stochastic behavior to the deterministic behavior. 
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ak = 0.2 X ^ 


ak = 0.9 X a] 
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FIG. A9: Effects of noise persist in a model of the SynExSlow circuit. Top row shows the deterministic ComK time 
series from the SynExSlow model (see Sec. A 2), while the next three rows show the stochastic ComK, MecA, and ComS time 
series for the same model. Although the SynExSlow model only exhibits a damped oscillatory regime at these parameters, 
not a standard oscillatory regime (see Sec. A 2), we define a heuristic boundary = 0.15/s below which oscillations are not 
appreciably damped within the first 24 hours (column 3), and above which they are (column 4). We see that, as in Fig. A6, 
noise prevents damping at large values of the control parameter, even at high molecule numbers (column 4). However, as in 
Fig. A8, noise does not induce repeated excitations at small values of the control parameter (column 2). We conclude that the 
former effect is more robust. 






































































